Introduction

Biological invasion refers to the process by which a living organism is transmitted from one place to another new environment through natural or artificial modes, causing economic loss or ecological disaster to the invasive biodiversity, agriculture, forestry, and human health (Mack 1997). Sicyos angulatus, also known as star cucumber, bur cucumber, wall bur cucumber and calabacilla (EPPO 2010; Osawa et al. 2016), is a native plant to the USA. As an ornamental plant, it has been imported into many European countries. Because of its greening effect of climbing growth and rapid morphing, S. angulatus has been widely used as a hedge plant and currently distributed in 36 states in the U.S.A. However, owing to its rapid expansion of encroachment ability, S. angulatus hinders the growth of other plants and even causes death. Delaware, Indiana, and Kenai define S. angulatus as a harmful weed to limit its spread. Besides, S. angulatus is also distributed in Japan (Kurokaw et al. 2009), Canada, Mexico, Korea (Kang et al. 2003) and the Netherlands (EPPO 2010). In Japan, S. angulatus was first discovered in Shizuoka Prefecture in 1952 (Kurokaw et al. 2009). The plant occurs in maize, peanut and soybean fields, and can also grow vigorously in non-cultivated land. In 2005, S. angulatus was identified as one of the most harmful weeds. From 2003 to 2013, S. angulatus was also found in China in Dalian, Liaoning and Hebei Provinces, including Beijing.

The preliminary research of S. angulatus have focused on its growth prevention techniques, external morphological characteristics, biological characteristics, hazard status, risk prediction and chemical prevention and control methods. However, S. angulatus seed germination characteristics and transcription profile related to germination and seedling growth, including associated metabolic pathways, have been rarely reported. In the present study, the transcriptome of S. angulatus was de novo assembled and differentially expressed genes (DEGs) in various growth stages of germinated seeds and seedlings were detected by transcriptome sequencing technology. The results obtained will provide evidence of invasive biology and a theoretical basis for efficient control strategies for S. angulatus invasion.

Materials and Methods

Plant material

The seeds of S. angulatus were germinated in an incubator at 25°C, during which water was added occasionally to retain moisture. After 7 days, three groups of germinating seeds without germination standard, with germination standard and cotyledon growth stage, respectively, were collected and stored in a freezer at −80°C after quick freezing in liquid nitrogen.

RNA extraction and sequencing

The total RNA was extracted by the phenol/chloroform method and detected by 1% agarose electrophoresis. A Kaiao K5500 spectrophotometer and Agilent 2100 RNA Nano 6000 Assay Kit were used for RNA quality assessment. After the total RNA samples were assessed, Oligo (dT) magnetic beads were used to enrich the mRNA, and fragment buffer was added to the obtained mRNA to generate short fragments that were used as templates to synthesize the first strand of cDNA with six base random primers. The second strand of the cDNA was further synthesized by adding buffer, DNA Polymerase I, dNTPs and RNaseH, purified by QIAQuick PCR kit. Next, cDNA was subjected to terminal repair by adding base A and sequencing adapter, and then selected by agarose gel electrophoresis. The final fragment was amplified by PCR, and the nine libraries obtained were sequenced by Illumina HiSeq™ 2000.

De novo assembly

First, the adapter and low-quality sequences were removed from the raw reads, including reads with N percentage of >5%. Second, the clean reads were used for de novo assembly and mapping to the transcriptome. Trinity (Grabherr et al. 2011) assembly program was employed for de novo assembly with default parameters. Third, the obtained unigenes from nine libraries were further spliced and assembled to obtain non-redundant unigenes. Finally, Bowtie2 (Langmead and Salzberg 2012) was employed to map the clean reads to the assembled transcripts sequence, and the ratio of the mapped reads was counted.

Functional annotation

TransDecoder software in the Trinity software package was used to identify the open reading frame (ORF) region of the assembled unigenes (Grabherr et al. 2011). The criteria for identifying ORF were as follows: ORF length >200, log-likelihood score >0 and if one ORF contains the other, the longest one is the output. To obtain comprehensive information on gene function, the predicted ORF was annotated by using Trinotate. BLAST and HMMER were also used to search the NCBI, Uniprot, evolutionary genealogy of genes: Clusters of Orthologous Groups of proteins (COGs), Non-supervised Orthologous Groups (eggNOG) and PFAM databases, respectively. SignalP , TMHMM (v. 2.0c) and RNAmmer (v. 1.2) were employed to predict the signal peptide and rRNA, respectively (Lagesen et al. 2007; Nielsen 2017). Blast2GO software was used to obtain the Gene Ontology (GO) terms of the unigenes (Götz et al. 2008).

DEGs analysis

We used Bowtie2 to map the clean reads back to the transcriptome with the default parameters. The number of mapped reads for each unigene was counted and then normalized into RPKM value (Reads Per kb per Million reads). DESeq2 was employed to identify DEGs among the samples, and control of false positives was achieved by multiple tests and correction based on Benjamini and Hochberg methods, with q < 0.05 and |log2_ratio|≥1 indicating DEG (Love et al. 2014).

Fig. 1: Length distribution of contigs (A) and unigenes (B)

GO and KEGG enrichment analysis

GO and KEGG enrichment analyses were based on the GO and KEGG terms of DEGs and all the unigenes. The significances were determined based on the hypergeometric test and multiple tests with q < 0.05. For each GO term, at least five genes were mapped.

Results

Sequencing and assembly quality statistics

To obtain a general overview of the S. angulatus transcriptome profile during seed germination, the total RNA samples were isolated from ungerminated seeds (NG), germinated seeds (G) and two cotyledon seeds (TL), respectively. Subsequently, cDNA libraries were prepared and three independent Illumina sequencing runs generated a total of 491,967,468 short sequence reads consisting of 61,495,933,500 nucleotides (nt) in total, with an average length of 125 bp for each short read. After strict quality filtering, 162 million, 166 million and 162 million reads from NG, G and TL samples, respectively, were used for de novo assembly (Table 1). We used Trinity assemblies for further analysis, which generated 185,654 contigs with N50 of 1143 bp and 127,874 unigenes with N50 of 807 bp (Table 2). The contig and unigene length distributions are shown in Fig. 1. After assembly, bowite2 was used to map the RNA-Seq reads to the contigs, and the mapping results showed that 412,857,414 (83.92) paired reads can be mapped to the contigs (Table S1), suggesting the high quality of assembly.

Functional annotation and classification of the S. angulatus transcriptome

Table 1: Reads statistics of different transcriptome samples

Samples	NG_1	NG_2	NG_3	G_1	G_2	G_3	TL_1	TL_2	TL_3
Raw Reads Number	74,599,058	68,813,124	70,359,372	66,266,126	72,200,576	73,878,140	68,450,544	68,380,202	67,788,690
Clean Reads Number	55,456,416	53,520,204	53,767,774	53,149,730	55,731,946	57,533,782	54,685,270	53,528,424	54,593,922
Reads Length (bp)	125	125	125	125	125	125	125	125	125
Q30 Bases Rate (%)	94.15	94.62	94.27	95.32	94.68	94.88	95.57	95.01	95.91

Table 2: Sequence statistics of contigs and Unigenes

Basic Stat	Trinity	Unigene
N50	1143	807
N90	302	267
Min	201	201
Max	15823	15823
Count	185654	127874
Mean	726.9678919	605.1526737

Fig. 2: Characteristics of homology search of S. angulatus unigenes. (A) Length distribution of ORFs; (B) Venn diagram for number of unigenes annotated by BLAST with an E-value threshold of 10^-5 against protein databases; (C) GO annotation of unigenes

To annotate the transcriptome of S. angulatus, Trinotate was used for unigene annotation. Based on Uniprot, NR and NT databases, 44,660 unigenes were annotated in total (Fig. 2A). About 71.61% of unigenes could be annotated with sequences of Cucumis spp., such as Cucumis melo (17,047) and Cucumis sativus (14,935). Next, TransDecoder was employed for ORF prediction. In total, 35,176 coding sequences were predicted and translated into protein sequences, of which >31.27% (10,998) had a length of >300 AA, with 611 longest unigenes having a length of >1000 AA (Fig. 2B). Based on sequence homology, these unigenes were categorized into 67 functional groups belonging to three main GO ontologies: “Molecular Function,” “Biological Process,” and “Cellular Component”. As shown in Fig. 2C, the GO results revealed that a high percentage of genes from categories, such as “cellular process,” “metabolic process,” “cell part,” “organelle,” “binding,” and “catalytic,” played essential roles in plant growth and development, with only a few genes related to “cell aggregation,” “symplast,” and “morphogen”. It is noteworthy that the above-annotated unigenes were involved in almost all the main functions of S. angulatus. In addition, the non-homologous unigenes functions were annotated using BLAST against the COGs database. The results revealed that 11,734 unigenes had COG classifications, and belong to 24 COG categories. Among these COG categories, “Translation, ribosomal structure, and biogenesis” was represented the largest group (1748, 14.90%), followed by “General function prediction only” (1397, 11.91%), “Posttranslational modification, protein turnover, chaperones” (1187, 10.12%) and “Replication, recombination and repair” (1157, 7.57%). The following categories denoted the smallest groups: “RNA processing and modification” (56; 0.48%), “cell motility” (10; 0.09%) and “Nuclear structure” (4, 0.03%) (Fig. 3). These results suggested that the de novo assembled unigenes had wide transcriptome coverage.

Expression analysis during seed germination

The RPKM was used to evaluate the expression level of unigenes in S. angulatus. A total of 127,860 unigenes were expressed during seed germination and the average correlation coefficient of the three replicates was >0.9 (Fig. S1, Table S2). In addition, the densities of the nine libraries were plotted, and similar density curves were obtained (Fig. 4A). The findings indicated that the genes of the nine libraries had a similar expression profile. Furthermore, 48,591, 21,153 and 54,268 unigenes showed differential expression between G and NG, G and TL and NG and TL, respectively. When mixed together, 66,664 unigenes were differentially expressed during seed

Fig. 3: COG annotation of non-homologous unigenes

Fig. 4: Differentially expressed unigenes (DEGs) in the samples. (A) Unigenes expression distribution in the nine samples; (B) Venn diagram for number of DEGs in each sample; (C) Hierarchical cluster analysis of DEGs

germination and 8919 unigenes were similar in pairwise comparison (Fig. 4B). Subsequently, the DEGs were subjected to hierarchical cluster analysis (Fig. 4C). The result showed that the DEGs could be divided into two groups, one was the ungerminated (NG), and the other was germinated (G and TL), indicating that the transcriptome of the ungerminated seeds may be considerably different from the germinated seeds.

Functional annotation of DEGs

To better understand the biological functions of DEGs, the role of germination responsive genes was investigated. First, Blast2GO was used to obtain the GO terms of each unigene. Subsequently, the proportions of the upregulated and downregulated unigenes were classified into three main functional categories (Fig. S2). Cell part, cellular process, metabolism, binding, and catalytic activities were the major responsive classes that had the maximum number of unigenes, among which the number of upregulated unigenes was higher than downregulated genes. This result suggested that metabolic processes were enhanced; the binding and catalytic activities were increased during seed germination. Subsequently, the enriched GO terms of DEGs were analyzed. For the DEGs in G and NG (Fig. 5), terms related to photosynthesis were enriched, such as “photosynthesis, light harvesting,” “photosynthesis, light reaction,” “thylakoid part,” and “photosystem II.” In addition, terms related to the germination process were also enriched, such as “response to hormone,” “shoot system development,” “post-embryonic development,” “meristem development,” and “phyllome development.” For the DEGs in G and TL (Fig. 6), the terms “tissue morphogenesis,” “seed germination,” “post-embryonic development,” “embryonic root morphogenesis,” “photosystem,” “photosystem I,” “photosystem II,” and “photosystem II oxygen evolving complex” were enriched. For the DEGs in NG and TL (Fig. 7), the terms “shoot system development,” “response to hormone,” “regulation of meristem growth,” “meristem maintenance,” “regulation of meristem development,” “post-embryonic development,” “photosynthesis,” and “photosystem” were enriched. These results suggested that genes related to post-embryonic development, meristem development, and photosynthesis were differently expressed in the germinated seeds (G and TL), contributing to seed germination and maintaining meristem development. Besides, these findings also implied that the metabolic pathways of the germinated seeds were enhanced and that the hormone signal transduction processes were increased during S. angulatus seed germination.

Genes related to plant hormones

In the present study, genes related to plant hormone signal transduction pathway exhibited different patterns of induction among the samples. In plants, auxin is a key regulator of development (Liu et al. 2007). With regard to the auxin influx carrier protein, 14 and 17 genes encoding this protein family were found to be upregulated in the germinated seeds G and TL, respectively (Fig. 8). Furthermore, genes encoding auxin-responsive protein were also upregulated in G and TL (19 and 20 genes, respectively). In addition, four genes encoding the auxin response factor (ARF) were upregulated in the germinated seeds, whereas the gene c138936_g2, which encodes the ARF, was downregulated in TL. It must be noted that ARF is the transcription factor that activates and represses the auxin response genes. These results indicated that genes related to auxin signal transduction play important roles in seed germination.

Fig. 5: GO enrichment analysis of the DEGs between G and NG

Fig. 6: GO enrichment analysis of the DEGs between G and TL

Fig. 7: GO enrichment analysis of the DEGs between NG and TL

Gibberellin (GA) is an important endogenous hormone that has diverse effects on plant growth and development, and its seed germination promoting effects have been proved in many plants. In this study, genes related to GA signal transduction were upregulated in the germinated seeds (Fig. 8), such as the gene encoding GA-INSENSITIVE DWARF1 (GID1, c126089_g1), four genes encoding DELLA proteins (c145425_g1, c101933_g1, c134844_g1 and c134844_g2) and downstream transcription factor phytochrome-interacting factor 3 (PIF). While PIF4 was only upregulated in TL, both PIF3 and PIF4 were downregulated in G. These results suggested that PIFs are important for seed germination, but their functions reduce after seeds germination.

Fig. 8: KEGG annotation of DEGs related to auxin and GA signal transduction

Discussion

In recent years, S. angulatus L. has obvious adverse effects on maize production and biodiversity in Beijing, Liaoning, and Hebei, China. The fruit of S. angulatus L. has numerous spines, which allow its attachment to animals resulting in natural dispersal, and the seeds of the plant are dispersed by water. To prevent further invasion and reproduction of S. angulatus, the seed germination process of it were studied using RNA-Seq. The De novo assembly was performed using 491,967,468 reads from ungerminated seeds, germinated seeds, and two-cotyledon seeds transcriptome assemblies with Trinity pipeline, and obtained a high-quality transcriptome of S. angulatus with 185,654 contigs with N50 of 1143 bp and 127,874 unigenes with N50 of 807 bp. In addition, the transcriptome assembled comprised numerous unigenes, and some of the unigenes contained more than one transcript. This may be because Trinity assembler generates high numbers of putative transcripts, including the reconstructing transcripts and alternatively spliced isoforms, and cannot completely ascertain the structural basis for the observed transcript variations, thus, leading to the generation of similar transcripts (Haas et al. 2013). In total, 44,660 unigenes were annotated in the Uniprot, NR, and NT databases. The majority of unigenes exhibited significant homology to cucurbitaceous plant sequences. This was coincident with the species of the S. angulatus what it was belong to Cucurbitaceae family, and the alignment, to a certain degree, could reflect the degree of relatedness to the reference species in the NCBI. Most of the unigenes were homologous to those of C. melo, which also belongs to the Cucurbitaceae family (Garcia-Mas et al. 2012), indicating S. angulatus has a closely relationship with C. melo. Within the GO and COG function classification, the unigenes could be categorized into 67 functional groups belonging to three main GO ontologies and 24 COG categories. These unigenes were involved in almost all the main functions, such as “cellular process,” “metabolic process,” “cell part,” “Translation, ribosomal structure, and biogenesis,” and “General function prediction only,” indicating that the assembled unigenes had a wide transcriptome coverage. These results are similar with previous reports that the unigenes related to above term were dominated in the de novo assemble transcriptome (Cao and Deng 2017; Gao et al. 2018). In summary, the high-quality assembled and annotated of the S. angulatus transcriptome of may contribute to the study of its invasive mechanism and corresponding control measures.

Seed germination is the first and most important step for S. angulatus invasion. It is important to examine the DEGs of different germination stages of S. angulatus, which could reveal the molecular mechanism of seed germination and help in effective prevention of its invasion. In the present study, a large number of unigenes were identified as DEGs in ungerminated and germinated seeds, indicating that the transcriptome of seeds undergoes substantial dynamic changes during the seed germination process. In fact, the functions of DEGs were noted to be consistent with the biological process. For instance, genes related to the GO term “photosynthesis” were enriched in the germinated seeds, suggesting that the germinated seeds had initiated photosynthesis which could produce energy to meet the needs of development (Galili et al. 2014). In addition, genes related to “response to hormone”, “shoot system development”, “post-embryonic development”, and “meristem development” was also enriched in the germinated seeds. As we known, plant hormones, such as GA, auxin, brassinosteroids, and oxylipins play important roles during seed germination (Rajjou et al. 2012; Resentini et al. 2015; Chahtane et al. 2018). In the present study, we found the DEGs related to auxin and GA signal transduction pathways were enriched in germinated seeds, when compared those in ungerminated seeds. In addition, during seed germination, the shoot meristem system development is an important symbol of seed germination. We found ARF3 (c136566_g1) which is involved in regulating shoot meristem initiation in rice (Nagasaki et al. 2007), was up-regulated in germinated seed. With regard to DEGs in G and TL, genes associated with “tissue morphogenesis” were enriched. These findings indicated that not only hormonal variations, but also changes in morphogenesis and maintenance occur during seed germination.

Conclusion

The present study is the first to determine S. angulatus transcriptome in ungerminated and germinated seeds. The high-quality transcriptome of S. angulatus allowed us to examine the key genes and related pathways. To investigate the genes related to S. angulatus germination, 66,664 DEGs were identified in ungerminated and germinated seeds, whose functions well reflected the corresponding transitions during seed germination. Moreover, analysis of the auxin and GA signal transduction pathways revealed some key genes that were vital for seed germination. These results contribute to the study of S. angulatus and could help in preventing its further invasion.

Acknowledgements

JG and JZ conceived and designed the experiments. HS and YL participated in the analysis and wrote the paper. JD, WX and BZ collect the experimental materials. All of the authors read and approved the manuscript. This work was funded by the China Agriculture Research System (CARS-02-25); the project of Study on the Integrated Control Technology of New Invasive Alien Plants of Sicyos angulatus in Hebei Province (16226503D).

Author Contributions

Conflicts of Interest

The authors declare no conflict of interest.

Ethics Approval

Ethics approval is not required for this research.

Data Availability

It is declared that data relevant to this article are available with the corresponding authors and will be made available on demand.

Ethics Approval

Ethics approval is not applicable for this research.

References

Cao Z, Z Deng (2017). De novo assembly, annotation and characterization of root transcriptomes of three caladium cultivars with a focus on necrotrophic pathogen resistance/defense-related genes. Intl J Mol Sci 18; Article 712

Chahtane H, TN Nogueira Fuller, PM Allard, L Marcourt, E Ferreira Queiroz, V Shanmugabalaji, J Falquet, JL Wolfender, L Lopez-Molina (2018). The plant pathogen Pseudomonas aeruginosa triggers a DELLA-dependent seed germination arrest in Arabidopsis. Elife 7; Article 1-34e37082

EPPO (2010). EPPO data sheet on invasive alien plants: Sicyos angulatus. EPPO Bull 40:401‒406

Galili G, T Avin-Wittenberg, R Angelovici, AR Fernie (2014). The role of photosynthesis and amino acid metabolism in the energy status during seed development. Front Plant Sci 5; Article 447

Gao J-F, Y Gao, J-H Qiu, Q-C Chang, Y Zhang, M Fang, C-R Wang (2018). De novo assembly and functional annotations of the transcriptome of Metorchis orientalis (Trematoda: Opisthorchiidae). Exp Parasitol 184:90‒96

Garcia-Mas, J, A Benjak, W Sanseverino, M Bourgeois, G Mir, VM González, E Hénaff, F Câmara, L Cozzuto, E Lowy, T Alioto, S Capella-Gutiérrez, J Blanca, J Cañizares, P Ziarsolo, D Gonzalez-Ibeas, L Rodríguez-Moreno, M Droege, L Du, M Alvarez-Tejado, B Lorente-Galdos, M Melé, LM Yang, YQ Weng, A Navarro, T Marques-Bonet, A Aranda Miguel F Nuez, B Picó, T Gabaldón, G Roma, R Guigó, Casacuberta, M Josep P Arús, P Puigdomènech (2012). The genome of melon (Cucumis melo L.). Proc Natl Acad Sci USA 109:11872–11877

Götz S, JM García-Gómez, J Terol, TD Williams, SH Nagaraj, MJ Nueda, M Robles, M Talón, J Dopazo, A Conesa (2008). High-throughput functional annotation and data mining with the Blast2GO suite. Nucl Acids Res 36:3420‒3435

Grabherr MG, BJ Haas, M Yassour, JZ Levin, DA Thompson, I Amit, X Adiconis, L Fan, R Raychowdhury, Q Zeng (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29:644‒652

Haas BJ, A Papanicolaou, M Yassour, M Grabherr, PD Blood, J Bowden, MB Couger, D Eccles, B Li, M Lieber, MD MacManes, M Ott, J Orvis, N Pochet, F Strozzi, N Weeks, R Westerman, T William, CN Dewey, R Henschel, RD LeDuc, N Friedman, A Regev (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc 8:1494‒1512

Kang JH, BS Jeon, SW Lee, ZR Choe, SI Shim (2003). Enhancement of seed germination by aging, cold-stratification, and light quality during desiccation in burcucumber (Sicyos Angulatus L.). Kor J Crop Sci 48:13‒16

Kurokaw S, H Kobayashi, T Senda (2009). Genetic diversity of Sicyos angulatus in central and north-eastern Japan by inter-simple sequence repeats analysis. Weed Res 49:365‒372

Lagesen K, P Hallin, E Rødland, H Stærfeldt, D Ussery (2007). RT: RNammer: Consistent annotation of rRNA genes in genomic sequences. Nucl Acids Res 35:3100‒3108

Langmead B, SL Salzberg (2012). Fast gapped-read alignment with Bowtie 2. Nat Meth 9:357–359

Liu PP, TA Montgomery, N Fahlgren, KD Kasschau, H Nonogaki, JC Carrington (2007). Repression of AUXIN RESPONSE FACTOR10 by microRNA160 is critical for seed germination and post-germination stages. Plant J 52:133‒146

Love MI, W Huber, S Anders (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550–570

Mack RN (1997). Plant invasions: Early and continuing expressions of global change. In: Past Future Rapid Environ Changes, pp:205‒216. NATO ASI Series, vol 47. Springer, Berlin, Germany

Nagasaki H, J Itoh, K Hayashi, K Hibara, N Satoh-Nagasawa, M Nosaka, M Mukouhata, M Ashikari, H Kitano, M Matsuoka (2007). The small interfering RNA production pathway is required for shoot meristem initiation in rice. Proc Natl Acad Sci USA 104:14867‒14871

Nielsen H (2017). Predicting secretory proteins with signalP. In: Protein Function Prediction, Vol. 10, pp:59-73. Human Press New York, USA

Osawa T, S Okawa, S Kurokawa, S Ando (2016). Generating an agricultural risk map based on limited ecological information: A case study using Sicyos angulatus. Ambio 45:895‒903

Rajjou L, M Duval, K Gallardo, J Catusse, J Bally, C Job, D Job (2012). Seed germination and vigor. Annu Rev Plant Biol 63:507‒533

Resentini F, A Felipo-Benavent, L Colombo, MA Blazquez, D Alabadi, S Masiero (2015). TCP14 and TCP15 mediate the promotion of seed germination by gibberellins in Arabidopsis thaliana. Mol Plant 8:482‒485